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STRESS ANALYSIS . AND , STRESS INTENSITY FACTORS FOR - FINITE 
GEOMETRY SOLIDS CONTAINING RECTANGULAR SURFACE CRACKS 

by John P* Gyekenyesi and Alexander Mendelson 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland , Ohio 44135 

SUMMARY 

The line method of analysis is applied to the Navi er -Cauchy equations 
of elastic equilibrium to calculate the displacement field in a finite 
geometry bar containing a variable depth rectangular surface crack under 
extensionally applied uniform loading. The application of this method 
to these equations leads to coupled sets of simultaneous ordinary differ- 
ential equations whose solutions are obtained along sets of lines in a 
discretized region. Using the obtained displacement field, normal stresses 
and the stress intensity factor variation along the crack periphery are 
calculated for different crack depth to bar thickness ratios. Crack 
opening displacements and stress intensity factors are also obtained for 
a through-thickness, center cracked bar with variable thickness ^ The 
reported results show a considerable potential for using this method in 
calculating stress intensity factors for commonly encountered surface 
crack geometries in finite solids - 

INTRODUCTION 

The main goal of fracture mechanics is the prediction of the load 
at which a structure weakened by a crack will fail. Knowledge of the 
stress and displacement distributions near the crack tip is of fundamental 
importance in evaluating this load at failure. During the early develop- 
ment of crack mechanics most of the effort was focused on through- 
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thickness cracks which could be characterized as two-dimensional. How- 
ever, part-through cracks are the most common type of crack defect found 
in actual service conditions (ref. 1). 

Because of its inherent three-dimensional character, only limited 
analytical work has been done in the past on surface crack problems. 

Early theoretical solutions , for surface flaw problems usually involved 
the discussion of semi-circular or semi-elliptical cracks in semi-infinite 
solids (refs. 2 to 8) . For this reason, results for finite geometry 
stress intensity factors are usually reported in terms of magnification 
factors applied to some convenient reference solution (ref. 9). 

These correction factors are then to account for certain inherent finite 
dimensions of the problem at hand. However, due to the formidable 
task involved in obtaining solutions for semi-infinite solids, only 
the front and back face magnification factors have been successfully 
calculated in the past, with little consideration given to the finite 
length or width normally encountered in engineering applications. In 
addition, considerable scatter exists in the reported results as obtained, 
by different investigators (ref. 9). 

Recently, approximate solutions of the finite geometry surface 
crack problem were obtained by the boundary integral equation method 
(ref. 10) and the finite element method (ref. 11). These solutions, 
however, have serious limitations since for acceptable accuracy in three- 
dimensional problems, the number of elements needed is extremely large, 
taxing even today computers. 

An alternate semi- analytical method suitable for the elastic solution 
of the surface crack problem is the line method of analysis. Successful 



application of this method to finite geometry solids containing cracks 
has been demonstrated by Gyekenyesl and Mendelson (ref. 12). Although 
the concept of the line method for solving partial differential- equations' 
is not. new (ref, 13), its application in the past has been limited in 
the analysis of solids and structures. 

It is the purpose of this report to present a simple and straight- 
forward approach to the elastic analysis of the finite geometry surface 
crack problem. Since the first treatment of this problem by Irwin in 
1962, considerable analytical effort has been expanded towards the 
development of an acceptable solution, but progress to this end has 
been slow and limited. In this report, the line method of analysis is 
extended to part-through rectangular cracks in bars of finite dimensions 
under normal tensile loadings. Results are. obtained for various crack 
depth to bar thickness ratios and effects of variable bar length and 
thickness on crack opening displacements are summarized. 

REDUCTION OF THE NAVIER-CAUCHY EQUATIONS TO SYSTEMS 
OF ORDINARY DIFFERENTIAL EQUATIONS 

Within the framework of linearized elasticity theory, the equations 
of elastic equilibrium in terms of displacements are 
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where the body forces are assumed to be zero and the dilatation is 


Su ^ ^ 

3x 3y 9z 


( 4 ) 


For a finite geometry solid with rectangular boundaries we construct 
three sets of parallel lines (fig. 1(a)). Each set of lines is parallel 
to one of the coordinate axes and thus perpendicular to the corresponding 
coordinate plane. An approximate solution of equation (1) can then 
be obtained by developing solutions, of ordinary differential equations 
along the x~directional lines. As seen in the figure, there are a total 
of £ = NY X NZ such lines were NY is the number of lines along the 
y-direction and NZ is the number of lines along the z-direction in a 
given plane, respectively. We define the displacements along these 
lines as u^, ...» The derivatives of the y-directional displace- 
ments on these lines with respect to y are defined as v' 

v’ I , and the derivatives of the z-directional displacements with 
respect to z are defined as * * ** These displace- 

ments and derivatives can then be regarded as functions of x only since 
they are variables on x-directional lines* When these definitions are 
used, the ordinary differential equation along a generic line ij (a 
double subscript is used here for simplicity of writing) in figure 1(b) 



may be written as 
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V* 

and 


dv 

dy 


w 


dw 

dz 


Similar differential equations are obtained along the other x--directional 
lines. Since each equation has the terms of the displacements on the 
surrounding lines, these equations constitute a system of ordinary differ- 
ential equations for the displacements u^, ^^2’* ’ ** 

The set of Z second order differential equations represented by 
(5) can be reduced to a set of 2^ first order differential equations by 
treating the derivatives of the u^s as an additional set of i unknowns, 
i.e. defining 


u 


£+1 


du^ 

~dx * ’^Jl+2 


du^ 

dx 


» etc. 


(7) 
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The resulting 2Z equations can now be written as a single first order 
matrix differential equation 

II = A^U + R(x) (8) 

where U and R are column matrices of 2Z elements each and is a 
2ix2Z matrix of the constant coefficients appearing in equations (5) 
and (7). 

In a similar manner, to solve equations (2) and (3) ordinary 
differential equations are constructed along the y- and z-directional 
lines respectively. These equations are also expressed in an analogous 
form to equations (8) ; they are 

(9) 

( 10 ) 

Equations (8) to (10) are linear first-order ordinary matrix 
differential equations. They are^ however, not independent, but are 
coupled through the vectors R, S and T whose components are given by 
equations similar to (6). The elements of the coefficient matrices A^, 

A^ and A^ are all constants, being functions of the coordinate increments 
and Poisson’s ratio only. 

Noting that a second-order differential equation can satisfy only 
a total of two boundary conditions and since three-dimensional elasticity 
problems have three boundary conditions . at every point of the bounding 


I7 ■ V + 

f - A H + T(z) 
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surface, some of the boundary data must be incorporated into the surface 
line differential equations. Hence i conditions of normal stress and 
displacement are enforced through the constants of the homogeneous solu- 
tions while shear stress boundary data must be Incorporated into the 
differential equations of the surface lines. The application of the 
specified shear conditions permits the use of central difference 
approximations when surface line differential equations are constructed. 
The details of constructing these equations are found in reference 14. 

SOLUTION OP THE SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS 

The systems of ordinary differential equations (8) to (10) can be 
solved by any of a number of standard techniques. The method used 
herein was basically the matrizant or Peano-Baker method of integration 
(ref . 15) . For equation (8) the solution can be written as 

U(x) = 

with similar solutions for equations (9) and (10). U(0) is the initial 

value vector, determined from the boundary conditions. The conversion of 
given boundary data into required initial values is discussed in more 

detail in reference 14. 

A. X 

The matrizant el is generally evaluated by its matrix series* 

For larger values of x, when convergence becomes slow, additive formulas 
may be used. In addition similarity transformations can be used to diago- 
nalize the matrix These various techniques for improving the accuracy 

are discussed in detail in reference 14. 
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Since equations (8) to (10) and their boundary conditions are highly 
coupled, it is generally impossible to directly evaluate their solutions. 
Thus, a successive approximation procedure must be employed where 
assumed values must be used initially for the required unknowns. The 
cyclic resubstitutlonof the obtained solutions into the coupling 
vectors and the boundary conditions will usually converge to the correct 
solution, depending mainly on the accuracy to which the required matrizant 
can be evaluated. 

Once the displacement field in the body has been calculated and 
the successive approximation procedure has converged, the normal stress 
distributions can be obtained directly by using the stress-dlsplac,ement 
equations. The shear stresses, however, can be evaluated only through 
finite difference approximations for the required displacement gradients. 

NUMERICAL RESULTS 

Figure 2 shows a finite geometry bar containing a traction free 
rectangular surface crack. Because of the symmetric geometry and loading, 
only one-fourth of the bar has to be discretized as shown in figure 2(b). 

The solution of this problem was obtained by using two different sets 
of lines along the coordinate axes so that the convergence of the finite 
difference approximations could be checked. In a given direction, uniform 
line spacing was used in all computations although this is not absolutely 
necessary. In general, an attempt was made to use a finer grid along 
the direction of largest variable change. The crack edge location with 
respect to the imposed grid was assumed to be halfway between nodes 
specifying normal stress and displacement boundary conditions, respectively. 
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The successive approximation procedure required for decoupling the 

three sets of ordinary differential equations was terminated when the 

difference between successively calculated displacements at every 

-3 

point was less than a preset value (10 )• As expected, the convergence 

rate of this successive approximation procedure was greatly dependent 
on the initial guess for the required unknowns in the coupling vectors 
and boundary conditions, For maximum computer efficiency, displacement 
data obtained from the use of coarse grids was interpolated to obtain 
improved starting values for the computations involving the final spacing 
of lines. The required initial quantities for the preliminary coarse 
grid calculations were taken to be zero in our work. All calculations 
were performed on a UNIVAC 1106 computer, using double precision 
arithmetic. 

Selected results of the dimensionless surface crack opening dis- 
placements are shown in figure 3. Note that the crack opening increases 
rapidly with crack depth for 0.167 0.834, slightly exceeding even 

the surface crack displacement of a through-thickness crack at ^ = 0.834. 
The plane strain solution for a finite width center cracked bar is also 
shown in figure 3 for reference. Final displacement values in this 
report were obtained from a set of 100, 140 and 140 x,y and z-directional 
differential equations, respectively. 

The maximum crack opening displacement variation with crack depth is 
shown in figure 4. For shallow cracks, this variation is almost a 
linear function of crack depth, reaching its maximum value near - 0.90. 
Interestingly, for “ - 1.0 or a through-thickness crack, the maximum 
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crack opening displacement decreases slightly from this critical value. 

In order to show the singularity of the stresses, the y-directional 
normal stress in the crack plane is plotted in figure 5 for ^ = 0,5, 

The results clearly Indicate the singular nature of along the crack 
periphery. Note that in figure 5(a), this stress is maximimi near the 
cracked surface of the bar and is minimum at the z = 0 plane. Figure 5(b), 
which shows the stress variation along the horizontal edge of the crack, 
indicates that is maximum at x = 0 and minimum at x = 2.0, 

The variation of crack opening displacement as a function of bar 
length for a fixed crack geometry is shown in figures 6 and 7. As 
can be noted from these figures, the finite length of the bar causes a 
considerable increase in crack displacements for values of L < 4c. 

Similar results were obtained in (ref. 12) for a through-thickness , 
center cracked bar. 

STRESS INTENSITY FACTOR 

It is customary in fracture mechanics to describe the plane elasticity 
crack opening displacement as a superposition of three basic deformation 
modes (ref, 16). Since the problem shown in figure 2 has geometric 
symmetry and is symmetrically loaded, only the opening mode of crack dis- 
placement is obtained. In terms of the stress intensity factor for 
the opening mode K^, the plane elasticity crack displacements near the 
crack tip are given by (ref. 16) 


I - 2(1 - V) 
ly=0 " G 


K, 



V 


plane strain 


( 12 ) 



11 


K, 


'y=0 (1 + v)G I 


_K 

2-tr 


plane stress 


(13) 


where R here Is the distance from the crack edge. Note that R » z - (t-a) 
along the horizontal edge and R * (W-c) - x along the vertical edge 
of the crack. Since three-dimensional problems are neither In a state 
of plane strain nor in a state of plane stress ^ the definition of a 
stress Intensity factor for these problems must be first established. 

Note that by definition the plane stress and plane strain stress intensity 
factors are equal while the displacements are approximately 12.5 percent 
different for v => 1/3. Since most published solutions for stress in- 
tensity factors are based on plane strain assumptions, equation (12) is 
selected to calculate the stress intensity factor. Rearranging this 
equation so that the dimensionless crack opening displacements can be. 
used leads to 


Ev 

0 c 
o 


C^Ki 


T 

C 


(14) 


where 

. 4(1 - Vh 

^ ^o \J 2ttc 

A plot of equation (14) as \/r/c 0 can then be used to calculate Kj.. 

Since the crack opening displacement is a function of both the thickness 
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and width variables, the previously defined stress intensity factor 
varies in both the x and z-directions. It should be noted that this 
description of is completely arbitrary and that it is questionable 
if it has any real significance in three-dimensional elasticity problems. 
However, values of K^. are presented here so that a comparison is 
possible with published plane strain solutions (ref. 16), 

The stress intensity factor variation along the crack periphery 
for different crack depths in a given bar is shown in figures 8 to 14. An 
important conclusion is Immediately obvious from the results shown. The 
point of maximum stress intensity factor for the rectangular surface 
crack shown in figure 2 is at the end of the semi-major axis or near 
the surface for all the crack geometries with ~ > 0 ^ 25 . Furthermore, 
the difference between the stress intensity factors at the surface and 
at the end of the semi-minor axis increases with crack depth, provided 
that the other dimensions of the problem remain fixed. These results, 
therefore, indicate that even though the back face approaches the deep 
edge of the crack, the point of maximum stress intensity factor does not 
change. It should also be noted that the minimum stress intensity factor 
is always at the comer of the crack and its change with crack depth is 
minimal, in contrast to the maximum stress intensity factor which changes 
about 400% over the crack depth range Investigated. 

Figure 15 shows the stress intensity factor variation along the 
z^direction for the through- thicknes s , center cracked bar. Note that 
is maximum near the surface and the results are symmetric about the 
center of the bar. Similar results were obtained in (ref. 12) for a 
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somewhat different geometry problem. 

A plot of the stress intensity factors at the semi-major and semi- 

minor axes of the crack as a function of crack depth is shown in 

Si’ ^ 

figure 16* These results indicate that for shallow cracks, ~ < 0.25, 
the stress intensity factor is no longer maxlmuni at the surface but 
reaches its greatest value at the semi-ininor axis. Also note that 
for deep cracks, the maxiitium stress intensity factor at the end of the 
semi-major axis exceeds that obtained for a through^ thickness crack. 

By comparing the results in this report for the rectangular sur- 
face crack problem to those of a through-thickness crack in an identical 
geometry solid, a part through crack correction factor, can be ob- 
tained for conveniently expressing the calculated maximum stress Intensity 
factor. In terms of this correction factor, let the surface crack results 
be denoted by 

K == C- (K ) (15) 

^lAX T 

where Kj is the maximum stress intensity factor for the through -thickness , 
center cracked bar. From figure 15 for the geometry in question, 

= 3.65 0 A plot of in the above equation as a function of 

crack depth is shown in figure 17. As expected, is considerably less 
than one for shallow cracks, but for deep cracks, > 0.7, the value of 
exceeds unity. 

Since the solution for an elliptical crack in an infinitely large 
body is readily available and convenient to apply, it is customary in. 
three-dimensional fracture mechanics to express results for finite geometry 
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problems in terms of this solution and suitable magnification factors. 

The crack opening displacement variation with bar length in figures 6 
and 7 can be used to calculate a finite length magnification factor 
for the selected crack geometry shown In figure 11, A plot of equation 
(14), as the bar length is increased from L = 1.25, yields a 7.7% in- 
crease in Kj for L = 1.75, when compared to the infinitely long case. 
Hence, = 1.077 which, by the way, differs considerably from the 
plane stress finite length correction factor given in (ref. 17) for 
the through-thickness, central crack problem. Other commonly applied 

magnification factors are M , M,., M, and M , where these factors account 

w r D p 

for finite width, front face, back face and plasticity effects, 
respectively. Values of these coefficients have been obtained previously, 
on a limited basis (ref. 9), but direct comparison of these numbers 
with results in this report is not practical. 

A plot of the surface crack opening displacement as a function of 
bar thickness is shown in figure 18 for a through-thickness, center 
cracked bar. Note that this displacement increases with increasing 
bar thickness. 

Figure 19 shows the stress intensity factor variation across the 
bar thickness for the same center cracked bar, Islda's plane strain 
solution (ref. 17), corrected for finite length and width, is also 
plotted in this figure. These results indicate that the maximum value 
of at the surface remains constant with variation in bar thickness , 
while at the center of the bar, approaches Isida’s plane strain 
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solution. However, the obtained minimum value of is about 17% higher 
than Isida’s result, which would indicate that the plane elasticity 
finite geometry magnification factors are not equal to those obtained 
from three-dimensional calculations. It should also be noted that using 
the calculated displacements in equation (13) rather than equation (12) 
would result in a 12.5% lower value of than that shown in figure 19. 
Furthermore, the maximum displacements in figure 18 differ less than 
10% as expected, while for t = 5.5 in figure 19, varies from 
3.0 to 3.65 or about 22%. 

Although the stress intensity factors for these problems could 
be determined with reasonable accuracy, the associated type of singu- 
larities are difficult to evaluate because values of the normal stresses 
are needed within a distance of 0.05c or less from the crack edge. With 
the equal spacing of lines used in these examples, the minimum node 
location for these problems is about 0.06c. For this range of crack 
edge distance R, the singularity of the stresses is not defined. 

CONCLUSIONS 

The line method of analysis presented affords a practical straight 
forward way for analysis of three-dimensional crack problems, at least 
for bodies with reasonably regular boundaries. Because parts of the 
solution are obtained as continuous functions along the lines chosen, 
relatively good accuracy can be obtained with coarse grids. Results of 
the analysis include the displacements and normal stresses at every 
node inside the body from which the stress intensity factor variations 
were easily calculated. In addition it should be noted that the common 
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semi-elliptical surface crack problem could also be analyzed by 
merely changing the boundary conditions at certain nodes in the crack 
plane. Introduction of plasticity into the analysis could also be 
accomplished by changing the coupling terms in equations (8) to (10) , 
Since these have to be determined by an iterative process in any case, 
it would seem possible to solve the elastoplastic problem by a simple 
extension of the present method. Whether this approach is practical 
requires further investigation. 
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(a) Three sets of lines parallel to x-, y-, and z-coordinates and 
perpendicular to corresponding coordinate planes. 





(b) Set of interior lines parallel to x-coordinate. 


Figure 1. - Sets of lines parallel to Cartesian coordinates. 
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L = 1.75 u = l/3 



Figure 3. - Surface crack opening displacement variation 
as a function of crack depth for a rectangular bar under 
uniform tension. 
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Figure 4. - Maximum crack opening displacement varia- 
tion with crack depth for a rectangular bar under uni- 
form tension. 
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(a) DIMENSIONLESS y-DIRECTlONAL NORMAL STRESS 
VARIATION ALONG BAR WIDTH. 



lb) DIMENSIONLESS y-DIRECTIONAL NORMAL STRESS 
VARIATION ACROSS BAR THICKNESS. 

Figure 5, - Dimensionless y-directional normal stress 
distribution in the crack plane for a bar under uni- 
form tension containing a rectangular surface crack. 
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Figures. - Variation of stress intensity factor Kj along the 
crack periphery for a bar under uniform tension contain- 
ing a rectangular surface crack. 
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figure 9. - Variation of stress intensity factor Kj along the 
crack periphery for a bar under uniform tension contain- 
ing a rectangular surface crack. 
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L*1.75 y = l/3 
W*2.0 a/t* 0.389 

c = 1. 0 a/2c » a 292 

t = 1.5 
a =0,584 
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Figure 10. - Variation of stress intensity factor Kj aiong the crack 
periphery for a bar under uniform tension containing a rectangu- 
lar surface crack. 





Figure 13. - Variation of stress intensity factor Kj along the crack 
periphery for a bar under uniform tension containing a rectangu- 
lar surface crack. 
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Figure 16, - Stress intensity factor variation with crack 
depth for a bar under uniform tension containing a 
rectangular surface crack. 
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Figure 17, - Surface crack correction factor variation 
with crack depth for a bar under uniform tension. 




Figure 18. - Surface crack opening displacement varia- 
tion as a function of bar thickness for a center 
cracked rectangular bar under uniform tension. 
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Figure 19, - Stress intensity factor variation as a function of bar 
thickness for a center cracked rectangular bar under uniform 
tension. 
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